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Abstract 



< 

■ We present a domain decomposition approach for the computation of the electro- 
, magnetic field within periodic structures. We use a Schwarz method with transpar- 
ent boundary conditions at the interfaces of the domains. Transparent boundary 
conditions are approximated by the perfectly matched layer method (PML). To 
cope with Wood anomalies appearing in periodic structures an adaptive strategy to 

\ determine optimal PML parameters is developed. 

We focus on the application to typical EUV lithography line masks. Light prop- 
agation within the multi-layer stack of the EUV mask is treated analytically. This 
CN \ results in a drastic reduction of the computational costs and allows for the simula- 

' tion of next generation lithography masks on a standard personal computer. 

o 
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1 Introduction 



The fabrication of semiconductor chips is essentially based on an optical pro- 
jection system. The pattern on a photolithography mask is transfered onto 

Email addresses: schaedle@zib.de, lin.zschiedrich@jcmwave.com, 
sven . burger® j cmwave . com, kloseOzib . de, frank . schmidtOzib . de (Frank 
Schmidt). 

^ Supported by the DFG Research Center Matheon "Mathematics for key tech- 
nologies" in Berlin. 



Preprint submitted to J. Comput. Phys. 



2 February 2008 



k = (fci, k2, k-ji) Illuminating Plane Wave 
Ail- 



Line 



X2 



X3 



/ 



Multi Layer Stack 



Substrate 
Period a 



Fig. 1. Layout of an EUV lithography hne mask. The structure is periodicahy re- 
peated in xi direction and invariant in X3 direction. The illuminating light is a plane 
wave with an arbitrary wave vector k = (ki, A;2, A^s). 



the chip by optical projection. State of the art photolithography tools are op- 
erated with light of a vacuum wavelength A ~ 193nm, [5]. Currently under 
development are tools that use extreme ultraviolet light (EUV) with a vacuum 
wavelength A ~ 13nm. These future systems will contain optical components 
including multi-layer structures which serve as mirrors. A typical section of an 
EUV lithography line mask is depicted in Figure 1. The line mask is invari- 
ant in X3 direction and periodic with period a in xi direction. The multi-layer 
stack may consist of more than 100 layers. The thickness of each layer is below 
a wavelength. The incident wave is twofold oblique - oblique with respect to 
the mask plane and oblique with respect to the multi-layer structure (conical 
incident). The polarization of the incident field is arbitrary. 

In Section 2 we introduce the mathematical setting of the arising scatter- 
ing problem and derive the radiating boundary condition in terms of Fourier 
modes. Further, we show that the exterior Dirichlet as well as Neumann bound- 
ary value problem is ill-posed in the presence of so called Wood anomalies, [18]. 
To deal with the large computational domains we propose a domain decompo- 
sition method (Section 5). The overall structure is split into sub-domains. The 
multi-layer sub-domain is treated semi-analytically (c.f. Section 2.2) whereas 
the other subdomains are discretized by the finite element method utilizing 
the PML method to approximate transparent boundary conditions. The PML 
method goes back to Berenger, [2]. Convergence of the method was proven 
in [15,16] and [14] for non-periodic problems. As is shown in Section 3 the 
PML method fails for periodic domains in the presence of Wood anomalies. 
As a remedy we propose in Section 3.1 a new automatic adaption of the layer 
size and the spatial discretization within the PML which leads to a quasi 
infinite layer thickness in the presence of Wood anomalies. In Section 4 we in- 
troduce a variational formulation to couple the PML to the interior problem. 
In contrast to Elschner et al. [11,10] the electromagnetic field E = {Ei, E2, E3) 
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is discretized with higher order Whitney elements for the {Ei, E2) component 
and Lagrange elements of the same order in the component. This allows 
for the accurate evaluation of Fourier coefficients needed for the coupling to a 
multi-layer stack as well for the computation of the far field coefficients. 

The domain decomposition approach for the wave equation was first studied 
for the scalar Helmholtz equation. Despres and Shaidurov proposed to balance 
the energy fiuxes across domain interfaces for Helmholtz problems, [9,26] . This 
idea was further expanded, [1,8,7,13,6,12]. Toselli used the PML method at 
the interfaces of the sub-domains, [27]. This idea is closely related to the ideas 
of multiple scattering, c.f. [17] with further references. In each sub-domain a 
simplified scattering problem is solved and the scattered field is added to the 
incoming field for the neighboring domains. We have recently presented an 
additive Schwarz algorithm for Helmholtz scattering problems with transpar- 
ent boundary conditions at the domain interfaces [20]. In this publication we 
used the DtN operator (Dirichlet to Neumann map). Since the definition of 
the DtN operator relies on the solvability of the exterior Dirichlet problem a 
DtN operator may not exist for periodic structures. Hence we avoid the usage 
of the DtN operator for the formulation of the domain decomposition method 
in this paper. 

The generalization to three dimensional geometries is straightforward. 



2 Scattering off periodic line masks 

The scattering off a periodic line mask is described by a Maxwell scattering 
problem, with Bloch-periodic boundary condition in Xi direction and transpar- 
ent boundary conditions in X2 direction. The dependency on the X3 component 
is eliminated. 

We consider electromagnetic scattering problems governed by the time-harmonic 
Maxwell's equations 



which may be derived from the Maxwell's equations when assuming a time 
dependence of the electric field as E(a;, t) = E (x) exp(— zcjt) with angular 
frequency uj. The dielectric tensor e and the permeability tensor /i are L°° 
functions of the spatial variable x = {xi,X2,X3). In addition we assume that 
the tensors e and /i do not depend on X3, that they are periodic functions in 
xi with period a, i.e. e{x + (a, 0, 0)) = e{x), ^{x + (a, 0, 0)) = fi{x), and that 
they are constant for X2 > X2,+ and X2 < X2- with X2,+ > X2-. For simplicity 



curl /i ^ {x) curl E (x) - ou'^e (x) E (f) = 0, 

dive (f) E (x) =0, 



(la) 
(lb) 
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assume that the dielectric and the permeabihty tensors are isotropic so they 
may be treated as scalar valued functions. Recall that any solution to (la) 
with 7^ also meets the divergence condition (lb). 

A scattering problem may be defined as follows: given an incoming electric 
field Einc satisfying the time-harmonic Maxwell's equations (la) for X2 > X2,+ 
and X2 < X2,-, compute the total electric field E, which satisfies (la) in M^, 
such that the scattered field Esc = E — Einc defined for X2 > X2,+ and X2 < X2~ 
meets the radiation condition given in Section 2.1. From a physical point of 
view, the scattered field has to be outward radiating, so it only transports 
energy towards infinity. 

It is possible to restrict the problem onto a two dimensional strip [0, a] x M 
provided that the incoming field is Bloch periodic in xi, [4] and depends 
harmonically on X3, i.e. 

Einc (xi + a, X2, xa) = Einc (xi, 0:2) e^'^^'^e^'^^"^ (2) 

where Einc is a periodic function in xi with period a. The important case of 
an incoming plane wave meets these restrictions. The total field E as well 
as the scattered field are then themselves Bloch periodic in xi and depend 
harmonically on X3. This can be seen using a symmetry argumentation where 
the unique solvability of the scattering problem is assumed. 

Below E, Einc and Egc denote the restriction of the respective field onto the 
strip [0, a] x M. Let us introduce the domains VL = [0,a] x [x2-,X2,+], ^+ = 
[0, a] X [x2,+ , 00] and accordingly With the definitions 



div geE = d.:,^eEi + d^^eE2 + iheE-^ 

the scattering problem splits into an interior domain problem 

curl3/i"^curl3E(xi,X2) -ci;^eE(xi,X2) =0 {xi,X2) eVt, 

E(0,X2) -E(a,X2)e''=^'^ =0, ^ ^ 

an upper exterior domain problem 

curl 3/x:^^ curl 3Esc,+ (xi,X2) - uj'^e+Esc,+ {xi,X2) =0 (xi,X2) G Vt+, 
E3c,+ (0,X2)-Esc,+ (a,a;2)e*'=^'^=0 

and a lower exterior problem on f2_ of similar type. 

Subproblems (3) and (4) are coupled by the following matching conditions on 
the boundary X2 = X2,+ 
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(E - (Esc,+ + Einc,+)) X f?+ = 0, (5a) 
(/i-^ curl 3E - (/i;^ curl 3Esc,+ + curl 3Einc,+)) x n+ = 0, (5b) 

where n+ = (0,-1,0)"^ denotes the unit normal vector. An analogous con- 
dition holds on the boundary X2 = X2,-, coupling the interior and the lower 
exterior problem. 

2.1 Radiation condition for homogeneous exterior domain problem 

The exterior domain problems lack a radiation condition. As both upper and 
lower exterior problem can be treated similarly, we consider the upper exterior 
domain problem only and drop the Without loss of generality and to simply 
matters, we assume that X2,+ = 0. 

Due to the periodicity of Egc exp(— ifciXi) the field has an expansion into 
Fourier modes 

Esc(xi, X2) = e+*^'^^'^ J2 en(a;2)e*"^"2-/^ (6) 

n&Z 

with the Fourier coefficients 

e„,(x2) = - /%-^^ifEsc(e,a;2)e-^«("2Va)^^ _ 
a Jo 

The field E„(xi, X2, X3) = en{x2) exp{i{n27r / a + ki)xi) expi^ik^x^) is a solution 
of Maxwell's equations (la) for X2 > 0. Hence inserting E„ in (la) yields 

E„ (xi, X2, X3) = 4^^e^("Wa+fci)xigifc2,„x2gifc3X3 + 

i{n2Tr/a+ki)xi -ik2,nX2„ik3X3 



with /c2,n = yko — (n27r/a + fci)^ — fef, where the branch cut of the square 
root is along the negative real axis and ko = uj^fjle. From this representation 
it is easily seen that the field can be decomposed into an incoming and an 
outgoing part. 

We have to distinguish three cases: 

(1) Re/c2,n > 0, Im/c2,n = [propagating mode) Both parts are propagat- 
ing plane waves with wave vectors {nl-Kja + ki, k2,n, k^) and (n27r/a + 
ki, —k2,n, ks) respectively. The second part transports energy in the —X2 
direction. We therefore require Cn- = 0. This corresponds to the well 
known Sommerfeld radiation condition. 
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(2) Re/c2,n = 0, Imfc2,n > {evanescent mode) The first part is evanescent in 
X2 direction while the second part increases exponentially. Therefore we 
again require Cn- = 0. 

(3) k2^n = {anomalous mode) In this case both parts are equal and constant 
in X2 direction. Energy is only transported in Xi and x^ directions. For 
the sake of a consistent notation we set e„^_ = 0. 

Hence the correct radiation boundary condition is e„,_ = for all n G Z, such 
that the Fourier coefficients of the scattered field are given by e„.sc = e„,+. 
The anomalous case is rare. For example for ki = and ^3 = it only occurs 
if a = 27r/ {k^n), hence the wavelength must be a multiple of the period a. 

2.1.1 Ill-posed exterior Dirichlet/ Neumann boundary value problems 

In our previous paper [20] the DtN operator was used to state the coupling 
between the different domains. However the DtN operator must not exist in 
the periodic setting - the exterior Dirichlet problem is ill-posed in the presence 
of anomalous modes. 

This may be seen by rewriting Maxwell's equations separated in Fourier modes. 
With kn = {ki + n27i/a, k2,n, ^3) the vectors e^^sc satisfy the algebraic relations 



kn ^ {kn ^ (^n,scj k^Cn^sc — 0, C^"^) 

en,sc-kn = 0. (7b) 

The first relation stems from Equation (la) and the second relation from 
Equation (lb). If j G Z corresponds to an anomalous mode, that is kj = {ki + 
j27i/a, 0, /ca), then ej^sc = (0, 1, 0) satisfies (7). Hence E^c — Cn.,sc Gxp{ikj ■ x) is a 
solution of the exterior domain problem with zero Dirichlet tangential bound- 
ary values. Therefore the Dirichlet boundary value problem is not uniquely 
solvable. Furthermore due to the divergence condition (7b) the vector ej^sc 
must be perpendicular to kj. Hence for boundary values with {di, 0, d^) -kj ^ 
the problem is not solvable at all. 

By an analogous argument one shows that the Neumann boundary value prob- 
lem is also ill-posed in the presence of anomalous modes. 

2.2 Scattering of an isotropic multi-layer stack - the Transfer Matrix method 

The Transfer Matrix method which according to [19] was developed by Schus- 
ter [25], will be reviewed here shortly. For more details, the reader is referred 
to [19,3]. 
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Suppose we are in the situation of Figure 1. Let us consider only the material 
stack with m finite layers, positioned at X2j, j = 0, . . . , m, with X2j < a;2j+i- 
For j = 1, . . . ,m the layer stack is given by the layer thicknesses X2j — X2j-i, 
and the material coefficients ej and fij. Additionally for the semi-infinite half- 
spaces we have eo, /io and e^+i, /^m+i- Since the Transfer Matrix algorithm 
is applied to each Fourier mode ki^n separately we drop the subindex n in 
this section. In each layer we define the local wave vectors kj = {ki, k2j, ks) 



and k j = (ki, —k2,j, k^) with k2j = \Juj'^ejfij — kf — k^, such that Re A;2,j > 

and ImA;2j > 0. For a given excitation ^ from above of the form Einc = 



semi-infinite half space a purely outgoing field is assumed, i.e Aq = 0. In 
the layers we have 6m unknowns - each A ot B has 3 components. In the 
lower semi-infinite domain the only unknowns are the three components of Bq 
of the purely outgoing field. In the upper semi-infinite domain there are six 
unknowns for the excitation and three for the reflected field. 
These unknown are determined by the following linear conditions arising from 
Maxwell's equations: there are 1 + 2 + 2m + 1 equations from the divergence 
condition. At the m + 1 boundaries of the layers there are 2(m + 1) matching 
conditions for the tangential components of the Dirichlet data and the same 
number of conditions from matching the Neumann data. 



Here Eq := Einc + Esc- The missing 4 conditions are the tangential components 
of the Dirichlet and Neumann data of the given incoming field. 
This yields a linear system of equations. To avoid large condition numbers due 
to the complex material tensors, in each layer ansatz functions with amplitude 
equal to 1 at the layer midpoint are used. 



3 Perfectly matched layer method 

In the previous section we discussed the homogenous exterior domain problem 
and derived transparent boundary conditions for each Fourier mode. Trans- 
forming back from Fourier space, this boundary condition would be non-local 
and somehow the anomalous case had to be treated separately. The perfectly 
matched layer method is an approximate transparent boundary condition, in- 
troducing only small reflections that are well under control. The reflections 

^ Here it is not assumed that the exciting field transports energy only in one direc- 
tion. 





Ej_i X n = Ej X n 
curlEj_i X n = fij curlEj x n 



at a; = Xj^i for j = 1, . . . , n -|- 1. 
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do not occur at the interface of the computational domain and the PML, but 
stem from the truncation of the PML. Another major advantage of the PML 
method is, that it fits in the finite-element framework, described shortly in 
Section 4.1 and thus does not introduce "full" blocks in the discretization. 

The PML method is based on a complex continuation of the scattered field. 
For 7 = (1 + ia), cr > 0, we define the complex continued field 

= ^ ^gi{n27ra+fci)xigifc„,27^2g«fc3a;3 j ^g^ 

With the definition 

curls^^E^ = {—d^^E^^s — ik3E^^2, iksE^^i — d^-^E^^^, dx^E^^2 — -dxi^-y^i) 

the field satisfies Maxwell's equations (4) with curls replaced by curl 3^.^. 
In the absence of anomalous modes E^ is evanescent for X2 ^ 00, 

|E^| < e-^^^C, 

with K = min„gz{Ini fcn,2, crRe fc„^2}- The idea is to restrict the complex con- 
tinued exterior domain problem to a truncated domain Qp = [0, a] x [0, p] 
and to impose a zero Dirichlet boundary condition at X2 = p- In case k is 
small or even 0, i.e. if we are "close" to an anomalous mode a special adaptive 
PML is used, where the thickness p is increased like 1/k, and the discretization 
points are distributed with an exponentially increasing mesh width guarantee- 
ing an effective discretization, c.f. Section 3.1. Thus the unbounded exterior 
problem (4) is replaced by the truncated exterior domain problem 

curl 3/i~^ curl 3E^,p(xi,a;2) - u;^£:+E^,p(xi, X2) =0 {xi,X2) G ^Ip, 

E^,p(0,X2)-E^,p(a,X2)e^'=^"=0, (9) 
^■y,p\x2=p X n =0. 

This modified truncated exterior problem is coupled to the interior problem 
using the modified matching conditions, c.f. (5) 

(E-(E^,, + Einc)) xn = 0, (10a) 
{e curl 3E — curl s^-yE^^p + curl sEinc)) x n = 0. (10b) 

Instead of homogenous Dirichlet boundary conditions at X2 = p in (9) one can 
equally well require that the Neumann data is homogenous, as in the PML 
the solution is oscillating and exponentially damped. 

3. 1 Automatic Adaption of PML 
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Algorithm 1. Adaptive PML method 



Require: e, a, /lint, 

Compute Np,„ and ,^max depending on hint and finite element order 
while (not converged) do 
eo = 0.0;ei = /iint;iV = l; 
while (-ln(e)/(^Arcr) < Kmin) do 

^N+i = + max{/iint, 2Tia^N/ (- ln(e))/A^p.w}. 
if (^7v+i > 1/e) then 
break 



N = N+1 
end if 
end while 

Compute solution u with PML discretization {^o, ■ ■ ■ , ^n} 

if \\u{;^N)\\<e\H-)\\ then 

converged 
else if > ^max then 

break 
else 



end if 
end while 



As discussed in Section 2.1 and Section 3 the PML method intrinsically fails 
in the presence of anomalous modes. For an anomalous mode the field behaves 
like exp{i{kiXi + /ca^s)) and hence a complex continuation in X2 direction has 
no effect on the decay property of the field. To obtain an effective transparent 
boundary condition we exploit the very specific behavior in X2 direction of 
the field and propose a mixed a priori and a posteriori refinement strategy 
of the perfectly matched layer method including the automatic adaption of 
the layer thickness p. The algorithm we propose is not restricted to the 2D 
periodic setting and was first published in [30] We therefore start from a simple 
generic model. As in our paper [23,29] a prismatoidal coordinate system in the 
exterior domain with a radial like coordinate ^ and an angular like variable rj 
is used. In the 2D periodic setting ^ is simply the X2 coordinate and rj = xi. 
Another example is a spherical coordinate system in 3D (r, (p, 9) with ^ = r 
and r} = (0, 9). To proceed we assume the following expansion of the field in 
the exterior domain 



with R,ek^{a) > and lmk^{a) > 0. Hence in ^ direction the field is a 
superposition of outgoing or evanescent plane waves. In the periodic setting 
such an expansion is explicitly given in (6). 



else 



■mm 





(11) 
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The complex continuation, ^1-^7^ with 7 = 1 + icr, gives 



Uy{v,0 II ~ / l|c(r/,«)||e 




(12) 



with K = crRe + Im k^. 

The PML method only effects the outgoing part with Re k^ strictly larger zero. 
Field contributions with a large Refc^ component are efficiently damped out. 
Furthermore evanescent field contributions are damped out independently of 
the complex continuation. For a proper approximation of the oscillatory and 
exponential behavior a discretization fine enough is needed to resolve the field. 
In contrast to that anomalous modes or "near anomalous" modes with fc^ ~ 
are neither evanescent nor damped out efficiently by the PML. Hence they 
enforce the usage of a large p but can be well approximated with a relatively 
coarse discretization in ^ due to their smoothness in ^. These requirements can 
only be satisfied by using an adaptive discretization. It is useful to think of the 
complex continuation as a high-frequency filter. With a growing distance ^ to 
the interior coupling boundary the higher frequency contributions are damped 
out so that the discretization can be coarsened. 

For a given threshold e selected according to the global accuracy requirements 
as described later we introduce the cut-off function 



With that at > each component in the expansion (12) with n > k,co,€{^') 
is damped out by a factor smaller than the threshold e, 

Assuming that this damping is sufficient we are allowed to select a discretiza- 
tion which must only approximate the lower frequency parts with k, < K,co,eiO 
for ^ > If we use a fixed number A^p.„ of discretization points per (gen- 
eralized) wavelength 27r/K we get the following formula for the a priori de- 
termination of the local mesh width h{^) = 2na / Kco,e{i) / N^,^. A good choice 
of A^p.w depends on the order of finite element used in ^—direction and need 
not to be adapted locally because the field depends smoothly in direction. 
Since Kco,e(0 ^ 00 ioi ^ ^ Q the local mesh width would be zero at = 0. As 
it is not reasonable to use a finer discretization in the exterior domain than 
in the interior domain we bound the local mesh width by the minimum mesh 
width /lint of the interior domain discretization on the coupling boundary, 

h{i) =max{/iint, 27ro-/Kco,,,,(0/^p.w}- 

The parameters e and A^p.w are also fixed accordingly to the interior domain 
discretization quality. The grid {'Co5G!'C25 . . . } is recursively constructed by 




ln(e)/e • 
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Fig. 2. Test problem for adaptive PML discretization. The lower material has 
an refractive index equal to ngub = 1-5, the upper material block consists of air 
{risup = 1.0). By Sneh's law the field is totally reflected for an incident angle equal 
to the critical angle "d^ = 180 • arcsin(1.0/1.5)/7r ~ 41.81. 




20 30 40 50 60 41 41.5 42 42.5 43 

angle of incidence angle of incidence 



Fig. 3. Left: Field energy error in the interior domain. The three lines (o, ' , +) 
corresponds to different refinement levels of the interior domain. Right: Zoom into 
left figure near critical angle. 

This way grows exponentially with n. To truncate the grid we assume that 
components in the expansion with k < Kmin can be neglected so that the 
grid {^0,^1, • • • ,^Ar} is determined by Kco,e,e(^Ar) < ft^min < Kco,eMN-i)- In the 
periodic setting there exists such a Kmin > in case no anomalous mode is 
present. 

As an a posteriori control we check if the field is indeed sufficiently damped 
out at ^jv, ||M(-,^Ar)|| < e||M(-)||. ^ Otherwise we recompute the solution with 
min/2 ^ . Since for an anomalous mode the field is not damped at all 
we restrict the maximum to < 7r//co/e- The pseudocode to the algorithm 
is given in Algorithm 1. 

To demonstrate the performance of the adaptive PML algorithm we compute 
the refiection of a plane wave at a material jump, c.f. Figure 2. We vary the 
angle of incidence from d = 20° to 'd = 60°. Further the incoming field is 
rotated along the 2:3 axis by an angle of 45°, so that the incidence is twofold 
oblique (conical). Hence the unit direction of the incoming field is equal to k = 



^ Here we assume homogenous Neumann boundary conditions for the truncation 
of the PML layer. If homogenous Dirichlet boundary conditions are chosen for the 
truncation of the PML layer, the sufficient damping of the Neumann data may be 
checked instead. 

^ This strategy proved useful in many experiment. However we consider to refine 
it. 
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20 30 40 50 60 20 30 40 50 60 

angle of incidence angle of incidence 



Fig. 4. Left: Thickness of the PML layer. At the critical angle the thickness is 
up to 10^ times larger than the diameter of the interior domain. Right: Number 
of discretization points used in the radial direction (X2). Although the needed 
thickness of the layer is huge the number of unknowns used in the PML layer remains 
moderate. 
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0.000205 


0.000820 
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0.000206 


0.000205 
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0.000051 
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Table 1 

Convergence of field energy at critical angle of incidence. The first column cor- 
responds to the interior mesh refinement step. The relative error of the electric 
field energy in the interior domain is given in column two, AE = |||Eex||^2 — 
||E/j||^2|/|||Eex||^2- In column three the relative error of the magnetic field energy 
AE' = I II curlEex||^2 — II curlEh||^2|/||| curlEex|| is given. For fixed PML thickness 
the solution converges as the interior mesh is refined. 

(cos 45° sin cos ^, sin 45° sin . The interior domain we used has as size of 
1.5 X 1 in wavelength scales. To measure the error we compute the field energy 
within the interior domain and compare it with the analytic value. In Figure 3 
the error is plotted for different refinement levels of the interior domain. The 
"+" line corresponds to the finest level. In Figure 4 the automatically adapted 
thickness of the PML is plotted (left) and the number of discretization points 
in ,^ direction (right). As expected a huge layer is used at the critical angle, 
whereas the total number of discretization points remains moderate. As can 
be seen in Figure 3 the maximum error appears at the critical angle. From 
that one may suspect a failure of the automatic PML adaption. But a closer 
analysis reveals that the chosen discretization in the PML layer is sufficient 
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as can be seen from Table 1. Here the thickness of the perfectly matched 
layer has been fixed and we further refined the interior domain. This way we 
observe convergence to the true solution but the convergence rate is halved 
at the critical angle. Hence the maximum error at the critical angle comes 
from an insufficient interior discretization. We conjecture that this is due to a 
dispersion effect. Since the wave is traveling along the xi direction it reenters 
the periodic domain leading to large "path length" . 



4 Variational form 



The coupled problem given by (9), (3) and (10) can be casted into a variational 
problem on the Sobolev space i^o.p ( curls, U fip) of -^^(curls) fields with 
generalized zero Dirichlet values at X2 = p- 

For a given test function $ G Hq p ( curl 3, U Qp) the following identity holds 
true, 



7 / $ ■ curl 3 .^/i ^ curl 3,-yE^ 

= 7/ curl 3_^$ ■ /i"^ curl 3_^E^ — / $ ■ curlsEsc x n, (13) 

J sip J X2=0 

where E^(a;i, a:2, X3) = Esc(a^i, 7x2, X3), c.f. (8). We first proof this identity for 
7 G M \ {0} . Using the non-euclidian coordinate change 

T~^ : (xi,X2,X3) (xi,7"V2,X3) 

and applying the transformation rules for differential forms, see [28], one gets 

/ curl 3/i^^ curl 3Esc = / curl 3/1"^ curl sE^. and (14a) 

/ curl 3$* yu"^ curl 3Esc= / curl 3$, ■ fi~^ curl ^E^ (14b) 

with 



fi^ = I J|J VJ * (15a) 

E^(xi, X2, 2:3) =J*Esc(a;i, 7x2,^3) (15b) 

^^{Xi, X2, x-i) = J^^{Xi, X2, x-i) = J*$*(xi, 7x2,^3). (15c) 

J = diag(l, 7, 1) is the constant Jacobian of T. Note that E,. are the pulled 
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back fields to and Egc in the sense of differential form calculus. 
We have 



7 / curla^^/i ^curl3^^E^= / ■ curl syU^ curl 3E* (16a) 
7/ curls^^^ • yU"^ curla^^E^ = / curls^* ■ yU^^ curlsE* (16b) 



which is verified by inserting (15) and using curl 3,^ = |J| '^JcurlaJ*. 
On the other hand integration by parts yields 



/ ■ curl ^ curl 3E* 

= / curl 3$* • curl 3E* - / ■ ('//"^ curl3E* x fi 



(17) 



and a respective equation for Egc, fi, and $* with the domain of integration 
Q^p. These together with equations in (14) give 

/ $7- f/i7^curl3E* X = / 1^- (/i"^curl3Esc X n) . (18) 

Jx2=0 ^ ' Jx2=0 ^ ' 

Using that E^ = J~*E* and using that the tangential components of are 
equal to one derives from (16) and (18) the desired identity (13) for real 
7. Since each term is a holomorphic function in 7 the identity (13) holds true 
for 7 G C \ {0}. 

The coupled problem given by (9), (3) and (10) in weak form is given by 

/ curl 3$ ■ /i"^ curl 3E — u?^ ■ eE + 

Jn 

7 y" curl3^^$ ■ /i"-^ curl3^^E.^ - ■ eE^ (19) 
= — $ ■ ( curl 3E — curl 3Esc) x n 

JX2=0 

Due to the Neumann coupling condition curlsE x n = curl3Esc x n + 
curlsEinc X n the boundary term is equal to /^.^^q ^ ' /^~^ curl3Einc x n. 
This is not yet the basis for a Galerkin ansatz in Hoplcurl^^QU Qp) as 
there is a jump of the Dirichlet data across the boundary X2 = 0, precisely 
E^ + Einc = E|a.2=o- Let n(Einc X n) G i^o.p ( curls, f2p) denote an exten- 
sion of a field with tangential Dirichlet data equal to Ejnc x n at X2 = to 
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Hqp ( curia, ^p) and add this to to obtain 

/ curl 3# ■ curl 3E — cu^^ ■ £:E + 
Jn 

7 / curl3_^$ ■ curl3_^(E^ + n(Einc X n)) - ■ e{E^ + n(Einc x n)) 
= - !• ■ cur^Einc X n) + 

7 J curl3_^$ ■ /i""^ curl3^^n(Einc x n) — u;^$ ■ £:n(Einc x n) . 

(20) 

This motivates the definition of the composed field u G Ho^p ( curls, ^ U ^p) 
by u\n = E and u|n = Egc + n(Einc x n) and of the following bilinear form: 



with 



a($, u) := an($b, u\n) + an^i^ln^, u|nj (21) 

af2($,u):= / curl3$ ■ /x"-^ curlsu - 0;^$ ■ eu , (22) 
Jn 

aQp($,u):=7^ curl3,^$ ■ curl3,^u - ■ eu . (23) 

6r($,*) := / . (24) 



With 



we end up with the variational problem: find u G i/o,p (curl3,f2 U Qp) such 
that for all $ G ifo.p ( curl 3,QU flp) 

a($,u) = aQ^($,n(Einc X n)) - 6($, curlsEi^c x n) . (25) 

Here we have avoided the definition of a DtN-operator. The total field is 
calculated as the solution of a coupled system (computational domain coupled 
to the PML), where the Dirichlet and Neumann data enter the equation on the 
"right-hand side" . If u is a solution of Maxwell's equations (1), the integration 
by parts identity can be rewritten using these bilinear forms as 

an(#, u) — 6($, — curlau X n) = (26) 

This formula will be useful to represent the Neumann data. Note that in (26) 
n is the "inward" normal with respect to Q. 



4.1 Finite element discretization 



To discretize the variational problem (25) we use vectorial finite elements 
on a triangular mesh in the interior domain and on a quadrilateral mesh in 
the PML. The three sub-meshes - lower PML mesh, interior domain mesh and 
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upper PML mesh - fit non-overlapping. In the PML we use a rectangular mesh 
[0, Xi,2, . • • ,a] X [x2,+ ,X2,+ + ^i, . . . ,X2,+ + ^n] where ^i, . . . ,^Ar are determined 
as described in Section 3.1. Since the Sobolev space ifo,p( curls, U fip) is 
isomorphic to Ho p{ curl 20? ^ U Qp) x Hq ^{Q U Qp) with the two dimensional 
curl operator curl 20(^1, ^2) = dx^U2 — 8x2^1 we use higher order Whitney 
elements to discretize the first and second component of the electric field 
and standard Lagrange elements for the third field component of the same 
order. This finite element space is also used for waveguide mode computations, 
c.f. [24] and the references therein. 

Bloch periodicity is enforced by a multiplication of basis functions associated 
with one of two corresponding periodic boundaries of the domain by the Bloch 
factor, c.f. [4]. An interior edge element function remains unchanged, c.f. Fig- 
ure. 5 (left). The support of a basis function associated with a periodic edge 
on the boundary consists of two triangles, c.f. Figure. 5 (right). The restric- 
tion of the basis function to the left triangle is defined as the standard shape 
function, whereas the shape-function on the right triangle is multiplied by the 
Bloch factor exp(i/cia). The construction of Bloch periodic Lagrange elements 
is similar. 





Fig. 5. First order edge elements on a simple grid. In the interior the tangential 
component is continuous across element boundaries. At the Bloch periodic boundary 
there is a phase shift. 



5 Domain Decomposition Method 



The idea for the Schwarz algorithm with transparent boundary conditions 
at the interfaces is to calculate the solution on every sub-domain separately 
using transparent boundary conditions and iteratively add the scattered field 
of each sub-domain to the incoming field for the neighboring sub-domains. 
The presentation here is restricted to the multiplicative Schwarz-algorithm. 
In its most general form the domain-decomposition algorithm is given in (27). 
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1^1 


^3 




^2 

















^2,p 





Fig. 6. Schematic sketch of the various domains and PMLs. Left: The computational 
domain Q is spht in three sub-domains. Right: The sub-domain ^2, with its three 
different PMLs. 



There E" denotes the nth iterate on sub-domain ^j^p^i is the PML domain 
to Qj at the interface to Qi. and by Qj^p we denote the PML domain to Qj at 
the interface to the exterior, c.f. Figure 6. 



set Ej- = for all j 
while not converged 
for all sub-domains j 
find Ej such that 

a,($,E,) = an^, ^($,n(Ei,, x n)) + ^ a^^. ^ n(E, x n)) 

i 

- 6r, ($, curlgEine X n) - curlgE^ x n) 

i 

V$ G ifo,p( curl 3, VLj U fij-p Ui ^-p^i) 



This algorithm requires the evaluation of Neumann data curl sEj x n along 
the boundary. Moreover at cross points, i.e at points given by Q.j^p^k^VLjpi ^ 
ioi i ^ k or fij.pjt fl fi^.p ^ 0, the "incoming" field is not a solution of Maxwell's 
equations, as it may have jumps. 

This difficulty maybe overcome by choosing sub-domains with a large overlap 
and by coupling the incoming field to the computational not at the boundary 
but at some additional artificial boundary in the interior. Here we do not pur- 
sue this strategy, but avoid cross-points, by dividing the computational domain 
horizontally in several sub-domains. Thus the sub-domains are arranged in a 
linear way. Each sub-domain has only two well separated boundaries neglect- 
ing the periodic boundary, and at most two neighboring domains. Inserting an 
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Computational domain 



Material 

1x2,0 



Outgoing wave ^ ^ Reflected wave 
Material 

~ Multi-layer stack ~ 



2^2,0 



Substrate 

Fig. 7. Decomposition of the problem into two infinite sub-domains. The scatter- 
ing problem is solved by the Finite Element Method in the upper domain and 
quasi-analytically in the lower domain. 

additional post-processing step, the Neumann-data can be evaluated weakly. 

set Ej = for all j 
while not converged 
for all sub- domains j 
find Ej such that 

aj($,Ej) = +an,^($,n(Einc x n)) + ^ a^^. ^ n(Ei x n)) 

^ (28) 
- 6rj (^, curlgEinc X n) - 2Z^r,„(^, curlgEi x n) 

i 

V$ e i?o,p( curl 3, Vtj U fij-p U, fij-p,i) 
for all sub domains j' 

^6rj_,($, curlgEj x n) -f6r, (^, curlgE^ x n) = -aQ^($,Ej) 

i 

In order to distinguish ioi k the contribution fepj^^l^; curlaEj x n) from, 
say 6r,j($, curlaEj x fi), it is required that there are no test functions that 
have a support in elements adjacent to Tj^k and F^y simultaneously. 



5.1 Schwarz algorithm for EUV 



For the special application - scattering off an EUV-line mask - one can make 
use of the "simple" geometry of the double layer stack that serves as a mirror 
employing the Transfer Matrix algorithm of Section 2.2. 

A simple situation is depicted in Figure 7. The upper domain contains the 
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mask line whereas the lower domain consists of the multi-layer stack and the 
lower substrate block. Instead of solving Maxwell's equations by the finite ele- 
ment method in the multi-layer stack, the incoming field is Fourier transformed 
and for each Fourier mode the Transfer Matrix algorithm is used to calculate 
the scattered field. This can even be simplified. If the tangential component 
of each Fourier-mode vector field is written as the linear combination of two 
linear independent polarizations, it is sufficient to compute the reflection coef- 
ficients of the multi-layer layer stack for each mode and each polarization only 
once. The number of Fourier mode ranges from = rimin to Umax- To determine 
these, we set kmax = 0.1 • 2n /h^ax, where h^ax is the maximum segment size 
of a finite element at the boundary. Then rimax is the greatest integer such 
that ki + rimax'^T^/O' < kmax, i-G and rimm is the greatest integer, such that 

kx f^min^'^ / ^ kmax- 



6 Numerical Examples 

An academic example 

The simple geometry of this example is depicted in Figure 8. It consists of 
three domains VLi and 1^2, each with a quadrilateral material inhomogeneity 
and VL^ a layer stack of four layers below VL2. The period is a = 1. The different 
shadings correspond to different materials as indicated. The permeability is 
equal to 1 everywhere. The permittivity is given by Si = 1.01, 62 = 1.52, 
£3 = 1.03, £4 = 1.54, £5 = 1.55, 56 = 1-06, £7 = 1-57, eg = 1-08. The semi- 
infinite top and lower strips, with refraction indices £9 = 1 and Eq = 1 are not 
shown. These are completely modeled by the PML method. Hence there are 
rather big jumps in the material coefficients at domain interfaces. 

The incoming field is a plane waves with wave vector ki^c = (1, —2, 1) and 
wave length, A = 0.84. The strength is Sinc = (1, 1, 1) x /cjnc/IKl, 1, 1) x /cindl- 

In the experiment the relative error is measured against the discrete solution 
obtained by solving the scattering problem on the whole domain. In solving the 
the scattering problem on the whole domain, the PML is chosen adaptively. 
These PML parameters are then fixed and used for all subdomains. Three 
cases are distinguished. 

(1) Schwarz algorithm with two domains (D2): One domain is Qi and the 
second domain is the union of ^2 and ^3. Thus the layers are discretized 
by finite elements. In Figure 8 this corresponds to the dark gray lines. 

(2) Schwarz algorithm with two domains (D2-EUV): One domain is Qi, the 
second domain is fl2- ^3 the layer stack is treated analytically and is like 
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Fig. 8. Material distribution (left) and magnitude of the electric field (middle) for 
a simple test problem. Convergence plot (right) for kiac = (1,— 2,1), A = 0.84 and 
different refinement levels. 



a boundary condition for 1^2. That is, if the subproblem on ^2 is solved we 
iterate internally between ^2 and ^3 and stop if the error is below 10~^ 
or after at most 100 iterations. In the domain decomposition algorithm 
only the number of iterations between Qi and Q2 is counted. In Figure 8 
this corresponds to the black lines. 
(3) Schwarz algorithm with three domains (D3): We are using a multiplicative 
Schwarz algorithm with three subdomains. Within one "iteration cycle", 
we first solve for Ei, then for E2 and finally for E3. In Figure 8 this 
corresponds to the hght gray lines. 



For these three cases the error versus the number of Schwarz iteration cycles 
is shown in Figure 8 (right). The experiment is performed for three different 
refinement levels, where "x" corresponds to the coarsest level with 5920 de- 
grees of freedom on the whole domain including the PML, the next finer level 
labeled with "*" is obtained by one uniform refinement of the initial grid and 
the finest level labeled with "□" by two uniform refinements of the initial grid. 
In case (D2-EUV) the error saturates at a level that clearly depends on the 
refinement of the interior grid. This behavior can be expected as the number 
of Fourier coefficients that are taken into account to couple the layer-stack an- 
alytically in the Schwarz iteration is inverse proportional to the mesh-width. 
In case (D2) and (D3) the error saturates at le — 14, which is close to ma- 
chine precision. This surprisingly good convergence behavior will be further 
analyzed in a subsequent paper. 
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# degrees of freedom 



Fig. 9. Left: Sketch of an EUV line mask. Right: Error versus the number of degrees 
of freedom in finite element mesh. The dashed error curve is obtained using the 
domain decomposition algorithm, decomposing the computational domain in two 
sub-domains (the line and the multi-layer stack) and treating the multi-layer stack 
separately. The solid error curve is obtained discretizing the whole computational 
domain. 

Real life EUV mask 



A schematic sketch of a more realistic EUV line mask is shown in Figure 9. 
There only three out of ten MoSi double layers are shown. The periodicity a 
is 40nm. The line made of silicon (Si) and the chromium absorber (Cr) have a 
width of 20nm and a height of 15nm. The first silicon layer's height is lOnm. 
Each molybdenum layer (Mo) has a height of 6nm and the subsequent silicon 
layers have a height of 8nm. The wavelength is 14nm. The permeability is 1.0 
everywhere. The permittivities are Emo = 1-69 + 0.016z, esi = 1-21 + 0.002z, 
Ecr = 1-43 + 0.24i and EAir = 1-0. 

Starting from a coarse mesh the grid is pre-refined to have at least 3, 4, 5, 6, 
7, 8, 9, 16 and 20 points per wavelength locally. The solution obtained with 
20 points per wavelength is taken as a reference solution to measure the error. 
We use a domain decomposition algorithm and decompose the mask into Qi 
(line, absorber, air) and Q2 (multilayer-stack). The multilayer-stack is treated 
analytically as described in Section 5.1. Additionally we are using a damping 
factor of 0.66 in the domain decomposition algorithm to speed up convergence. 
The PML is chosen adaptively as described in Section 3.1. 

Figure 9 shows the error versus the number of degrees of freedom in the fi- 
nite element grid including the PML. To obtain the solid line, the multi-layer 
stack is discretized using finite elements. Clearly, if the multi-layer stack is 
not discretized, but treated analytically and coupled to Qi in the domain- 
decomposition algorithm, the number of degrees of freedom is reduced drasti- 
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cally. The above calculations where performed on an AMD Opeteron PC with 
16GB of RAM. The arising linear system problems are solved with the sparse 
LU method PARDISO, [21,22]. 

This reduction of the number of degrees of freedom due to the domain de- 
composition approach, allows to compute realistic masks on standard 32-bit 
computers. 
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